Ground-State Properties of Quantum Many-Body Systems: Entangled-Plaquette 

States and Variational Monte Carlo 
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We propose a new ansatz for the ground-state wave function of quantum many-body systems on a 
lattice. The key idea is to cover the lattice with plaquettes and obtain a state whose configurational 
weights can be optimized by means of a Variational Monte Carlo algorithm. Such a scheme applies 
to any dimension, without any "sign" instability. We show results for various two dimensional spin 
models (including frustrated ones). A detailed comparison with available exact results, as well as 
with variational methods based on different ansatzs is offered. In particular, our numerical estimates 
are in quite good agreement with exact ones for unfrustrated systems, and compare favorably to 
other methods for frustrated ones. 

PACS numbers: 02.70.Ss, 05.50+q 



I. INTRODUCTION 

The study of the ground-state (GS) properties of quan- 
tum many-body systems is one of the most challenging 
tasks of theoretical physics. Exact results can be ob- 
tained numerically only for systems of relatively small 
size (i.e., few particles); this limitation is particularly 
severe, e.g., when studying phase transitions, wherein 
the emergence of long-range order can only be estab- 
lished by carrying out an extrapolation of the estimates 
to the thermodynamic limit. Two main numerical tech- 
niques, namely Density Matrix Renormalization Group 
(DMRG)i and Quantum Monte Carlo (QMC)2 have 
successfully been employed to investigate quantum spin 
models on a lattice. These computational approaches, 
however, find their optimal applicability under different 
specific constraints. DMRG yields extremely accurate re- 
sults in one dimension (ID) even for very large systems, 
but fails in describing the properties of the quantum GS 
in higher dimension due to the unfavorable scaling with 
the system size of the computational resources needed^ 
QMC, on the other hand, is the method of choice for 
quantum systems obeying Bose statistics in any dimen- 
sion, but suffers from the notorious "sign problem" for 
Fermi systems. 

Generalizations of the variational family of Matrix- 
Product States (MPS) underlying the DMRG have re- 
cently been investigated to go beyond the discussed 
limitations of DMRG itself and QMC. Specifically the 
most natural extension of MPS is given by Projected- 
Entanglcd Pair States (PEPS)^ which efficiently approx- 
imate ground states of local Hamiltonians 5 - and have 
been used to simulate, otherwise intractable, 2D quan- 
tum systems<&£ Other variational families of states have 
also been proposed and tested in 2D^ ' 10 i 11 : 12 i 13 ' 14 De- 
spite the promising results obtained so far, it seems very 
hard to use those methods in 3D [or even for systems with 
periodic boundary conditions (PBC)] due to the unfavor- 
able scaling of the required computer resources. 

A new possibility has recently emerged to combine the 



main advantages of DMRG and Monte Carlo in order to 
build new algorithms! 15 ' 16 Some of the m 15 i 17 can be used 
in 2D and it seems that the one based on String-Bond 
States^ may be used even for some 3D systems^ 

In this paper we introduce a new numerical technique 
that combines the strengths of QMC with an extension 
of PEPS to simulate lattices models, overcoming some 
of their limitations. We also test this technique with 
non-trivial models, and compare it with other techniques, 
including PEPS. Specifically, we propose a new class of 
states called Entangled-Plaquette States (EPS) which al- 
low an accurate characterization of the GS of quantum 
spin systems by means of a simple Variational Monte 
Carlo (VMC) algorithm (see Ref. [lj^ for a general re- 
view) . 

The basic idea underlying EPS can be schematically 
described as follows: Assume to cut a lattice in sev- 
eral sub-blocks (e.g. small PEPS) and extract for each 
sub-block the GS wave function. The wave function of 
the original system, expressed as the product of the sub- 
block wave functions yields reasonable (up to corrections 
which scale with the sub-block boundaries) estimates of 
GS energy and short-range (of the order of the sub-block 
size) correlations. These estimates dramatically improve 
if the sub-block size is increased and, more importantly, 
if overlapping sub-blocks (i.e., entangled plaquettes) are 
employed. The latter is the crucial point which allows, 
accounting for its correlated nature, a description of the 
quantum GS much more accurate than that obtainable 
with a simple non-overlapping-plaquette product state 
(i.e., in a mean-field fashion)^ A GS wave function 
whose weights are the product of variational parameters 
in one to one correspondence to the spin configuration 
of each entangled plaquette, appears, consequently, the 
natural choice for a variational ansatz. Therefore, our nu- 
merical approach is based on a family of states, namely 
EPS, which share many analogies with PEPS, and takes 
advantage of Monte Carlo sampling to estimate physi- 
cal observablcs of interest. It can be applied to systems 
of any spatial dimensionality and, as a pure variational 
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method (i.e., not involving imaginary time projection), 
is sign problem free. 

We test our numerical protocol on a variety of quantum 
spin models on a square lattice comprising as many as 
N = Lx L = 400 sites. Our energy estimates are in excel- 
lent agreement with those ( "exact" in practice) obtained 
by QMC for a system of lattice hard core bosons. In the 
presence of nearest-neighbor repulsion at the Heiscnbcrg 
point, we find an extrapolated (to infinite lattice size) 
value of the energy per site which differs from the QMC 
result^ by less than 2xl0~ 3 , and is more accurate than 
VMC estimates obtained with a Jastrow wave function^ 
In the case of a frustrated antiferromagnet (i.e., the so- 
called J\ — J2 model) , for which a sign problem exists in 
QMC, our energy estimates (whose error relative to the 
exact ones is less than 1.5 x 10 ~ 2 ) compare favorably with 
those obtained with PEPS, or fixed-node Green Function 
Monte Carlo (GFMC)3 



II. METHODOLOGY 

Consider a collection of TV spins | arranged on a L x L 
square lattice with N = L 2 (without loss of generality, 
here and in the following we will refer to this specific 
case). Provided a trial state \ip) = J2 n ^( n )\ n ) where 
n) = \ni, 122, . ■ . , njv) and rii = ±1 V i = 1, ... ,7V, the 
energy expectation value on the given state is: 
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E(n) 



P(n) = 
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and W(n) = W*(n) (real weights are assumed for sim- 
plicity). According to the variational principle, (E) is an 
upper bound of the GS energy which can be evaluated 
by minimizing Eq. [T]with respect to the weights. At this 
point we have to make an ansatz for the wave function: 
imagine to cover the lattice with plaquettes (say one of 
dimension l± x l 2 for each site) and assign a coefficient 
Cp p to all the possible 2 llXl2 spin configurations of any 
single plaquette. Given a global spin configuration |n), 
its weight can be expressed as follows: 
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W(n) = <n|V) = II °p 
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where Cp p depends only on the spin state of the np 
sites belonging to the Pth plaquette. With this choice 
the analytic expression of the derivative of Eq. [T] with 
respect to C" p is: 
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(4) 

The multidimensional summations in Eq. [1] and Q] can 
exactly be evaluated only for small N. For large sys- 
tems (N > 30) , one has to employ the Monte Carlo 
method. Specifically, the energy, as well as its deriva- 
tives, can be estimated from the same sample. More- 
over, the only quantity depending on the plaquette co- 
efficient with respect to which the derivative is taken is 
Dp(n P ) = pp^y gpff . By using Eq. [3] it turns out 

that Dp(np) is simply equal to 1/C" P (easy and fast 
to compute). Similarly to other works / 5 ' 16 , the steps of 
the basic variational algorithm adopted to calculate the 
GS energy are: i) Start from a randomly chosen initial 
configuration; ii) generate a large set of new configura- 
tions by flipping one or more spins via the Metropolis 
algorithm^ iii) evaluate the energy and its gradient vec- 
tor; iv) update all the Cp r s of a small step against the 
gradient direction; v) iterate from ii) until convergence 
of the energy is reached. It is worth mentioning that ex- 
pectation values of physical observables other than the 
energy can be evaluated according to Eq. [T] and [2] when 
H is replaced by a generic operator O. 

For a single-spin flip the acceptance probability is given 
by: 



A 



(5) 



where the flip is proposed for the jth spin and the prod- 
ucts run over all the plaquettes which include such a spin. 

In a typical calculation, we start with 2x2 plaquettes; 
once the energy has converged, the size of the plaquettes 
is increased to improve the estimate. The number of 
coefficients which need be stored in memory for a spin- 
1/2 system is N x 2 llxl2 , and can be reduced taking 
into account problem-dependent symmetries. For each 
optimization step, a few thousands updates are necessary 
to get a rough estimate of the energy and efficiently move 
the coefficients along the gradient. At the later stages of 
the simulation, to reach the optimal energy value, an 
important role is played by the gradient step which has 
to be carefully tuned. An example of how the error of 
the GS energy relative to the GFMC result decreases as 
a function of the plaquettes size is illustrated in Fig. [TJ 

The relative error is already small (less than 1%) for 
2x2 plaquettes and reaches a value < 10~ 3 when 4x4 
plaquettes are used. Numerical data refer to a system 
of lattice bosons (at half filling) which interact via an 
infinite on-site (hard core) repulsion on a 10 x 10 lat- 
tice with PBC. Since the total magnetization along z is 
a good quantum number, we performed the calculation 
in the canonical ensemble (i.e., in the S z = sector). 
Consequently we chose to update the configuration by 
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FIG. 1: (color online) Dependence on the plaquettes size of 
the error in the GS energy (computed with the method illus- 
trated in this work) relative to the GFMC result for a system 
of hard core bosons at half filling on a 10 x 10 square lattice. 
PBC are assumed. The dashed line is only a guide to the eye. 



flipping pairs of spins j and k for which S? = —S%. The 
expression of the acceptance probability for the pair up- 
date is a straightforward generalization of Eq. [5] 



III. RESULTS 

Estimates of the GS energy per site of a system of 
lattice hard core bosons, for three different lattice sizes, 
are presented and compared to GFMC results in Tab. 
HI Although the GFMC method is not purely variational 
(i.e., the GS wave function is projected out from a trial 
state via imaginary time evolution) our numerical data 
are in excellent agreement with GFMC ones, even for the 
largest system considered in this work (L = 20). 



TABLE I: GS energy per site (in units of the nearest-neighbor 
hopping integral i) for a system of hard core bosons onaLxi 
lattice with PBC. For L = 20, 4x3 plaquettes have been 
employed, 4x4 ones in the remaining cases. GFMC estimates 
have been also computed and are shown for comparison. 



L 


This work 


GFMC 


8 


-1.1004(4) 


-1.1008(3) 


10 


-1.0988(4) 


-1.0998(1) 


20 


-1.0952(7) 


-1.0967(4) 



As a second case study we investigate the J\ — J<i quan- 
tum spin Hamiltonian: 



H 



J\ Si ■ Sj + J2 S; 



(6) 
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where the first (second) summation runs over nearest- 
(next-nearest-) neighbor sites. For Ji = 1 and J 2 = 
the above Hamiltonian describes the antifcrromag- 
netic Hciscnbcrg model, which is isomorphic to a system 
of hard core bosons with nearest-neighbor repulsion of 



strength V=2t. This model has been extensively studied 
numerically in the past (see for instance Refs. [U and 
22l ). GS energies per site for various lattice sizes are pre- 
sented in Tab |TT] and compared to exact calculations^ 
and Stocastic Series Expansion (SSE) data^ 



TABLE II: GS energy per site of the antiferromagnetic 
Heisenberg model on a L x L square lattice with PBC. SSE 
estimates and exact results are also shown for comparison. 
Numerical values are given in units of Ji. 



L 


This work 


SSE a 


Exact 6 


4 


-0.7016(1) 


-0.701777(7) 


-0.7018 


6 


-0.6785(2) 


-0.678873(4) 


-0.6789 


8 


-0.6724(3) 


-0.673487(4) 




10 


-0.6699(3) 


-0.671549(4) 





"Ref. [10 
6 Ref. [51 



The energy per site is a monotonic increasing function 
of the lattice size. The simple formula^ 



E(N) =E OQ -f3 



AT 5 



(7) 



where (3 = 1.4377^i gives an extrapolated energy per 
site Eoo = —0.6683(3), which is in good agreement with 
the most accurate QMC results^ E^ E = 0.669437(5) 
and lower than other extrapolated values obtained by 
purely variational methods. For example, on the basis 
of a Jastrow wave function, Trivedi and Ceperley found 
E J J T ~ -0.6590.^ 

Calculations carried on with open boundary conditions 
(OBC) yield GS energies lower than those obtained with 
the general PEPS, or SBS method. For example, we get 
E = -0.6258(1) for L = 10 while the PEPS result is 
e peps = _Q.62515 and the SBS one E SBS 0.6225. 

Estimates of the spin-spin correlation function com- 
puted at the maximum distance on the lattice accord- 
ing to the formula Corr(L/2, L/2) = (S r • S r <) where 
r — r' = (L/2, L/2) arc shown in Tab lllll We find an ex- 
trapolated value of the staggered magnetization defined 
by M 2 (L) = Corr{L/2 7 L/2) = + b/L of 0.324(1) 
which is in reasonable agreement with M~~ = 0.3070(3) 



TABLE III: Spin-spin correlation function (computed at the 
maximum distance on the lattice) of the antiferromagnetic 
Heisenberg model on a L x L square lattice with PBC. SSE 
estimates (adjusted for different factors in the definition) are 
also shown for comparison. Numerical values are given in 
units of Ji. 



L 


This work 


SSE" 


4 


0.1807(4) 


0.17962(1) 


6 


0.1550(4) 


0.152568(9) 


8 


0.1432(4) 


0.13760(1) 


10 


0.1352(5) 


0.12855(2) 



"Ref. Ul 
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TABLE IV: GS energy per site (in units of Ji) for the J\ — J2 
model on a square lattice (with PBC) comprising 36 sites. 
Exact results are also shown for comparison. 



J2/J1 


This work 


Exact 11 


0.0 


-0.6785(2) 


-0.6789 


0.1 


-0.6377(1) 


-0.6381 


0.2 


-0.5985(1) 


-0.5990 


0.3 


-0.5616(1) 


-0.5625 


0.4 


-0.5277(1) 


-0.5297 


0.5 


-0.4985(2) 


-0.5038 


0.6 


-0.4860(2) 


-0.4932 


0.7 


-0.5255(1) 


-0.5300 


0.8 


-0.5843(1) 


-0.5865 


0.9 


-0.6453(1) 


-0.6491 


1.0 


-0.7091(1) 


-0.7144 


a Ref. [H 



TABLE V: GS energy per site (in units of Ji) for the Ji — J2 
model on a square lattice (with OBC) comprising 64 sites. 
PEPS results are also shown for comparison. 



J2/J1 


This work 


PEPS" 


0.0 


-0 61 567CS1 


-0 61 506 


0.1 


-0.5865(1) 


-0.5845 


0.2 


-0.5571(1) 


-0.5555 


0.3 


-0.5290(2) 


-0.5274 


0.4 


-0.5028(1) 


-0.5016 


0.5 


-0.4781(1) 


-0.4779 


0.6 


-0.4523(2) 


-0.4508 


0.7 


-0.4553(1) 


-0.4541 


0.8 


-0.4928(2) 


-0.4906 


0.9 


-0.5371(1) 


-0.5344 


1.0 


-0.5800(1) 


-0.5792 


"data provided by V. Murg 



reported in Ref [2l| and with other estimates^ It has to 
be mentioned that the discrepancy between our result 
and the SSE one might be due to a violation of the so- 
called "area law" i2£ 

By including the next-nearest-neighbor interaction 
( J2 > 0) in the J\ — J<x Hamiltonian the system becomes 
frustrated and is believed to undergo a phase transition 
for J2 — 0.6. This model cannot be simulated by QMC 
due to the sign problem (arising in turn from the under- 
lying Fermi statistics). Our method, instead, is purely 
variational and can be applied without the occurrence of 
any sign instability. With the aim of comparing our es- 
timates with exact results (available only for N < 36), 
we computed the GS energy, as a function of J2/J1, of 
the J\ — J2 model on a 6 x 6 square lattice. Numerical 
values are shown in Tab. IIVI the error of our estimates 
relative to the exact result is shown in Fig. [2j This 
quantity never exceeds 1.5xl0 -2 , moreover it has to be 
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FIG. 2: (color online) Error (relative to the exact result) of 
the GS energy of the Ji — J2 model computed with the method 
illustrated in this work. The square lattice comprises N = 36 
spins; PBC are assumed. The dashed line is only a guide to 
the eye. 



mentioned that the GS energies computed by means of 
EPS compare favorably to those obtained by GFMC with 
the fixed-node approximation^ and, except in a narrow 
region of J2/J1 ~ 0.5, where the the agreement is how- 
ever remarkable, to variational results obtained with a 
BCS-type ansatz, 2 ^ 

Also in this case, for OBC, we test our scheme against 
the PEPS one. The estimated energies are lower than 
those computed with PEPS even for a 8 x 8 lattices (see 
Tab. |V]) . where the PEPS approach performs at its best. 



IV. CONCLUSIONS AND OUTLOOK 

In this work we have presented Entangled-Plaquette 
States: an ansatz for the GS wave function of quantum 
many-body systems which allows the accurate estimate 
of physical observables by means of Variational Monte 
Carlo. Our approach not only gives accurate results for 
unfrustrated systems (where other methods are in prin- 
ciple "exact" ) but, most importantly, applied to systems 
for which QMC simulations suffer from a sign problem, 
yields estimates whose accuracy seems to be not obtain- 
able with different techniques. Therefore our method 
appears as a very promising avenue to investigate the 
effects of Fermi statistics (e.g., frustration). The exten- 
sion of our computational approach to 3D can be easily 
implemented (i.e., taking cubic plaquettes) and several 
possible improvements (e.g., cover the lattice with pla- 
quettes of different shapes and sizes) , are being currently 
investigated. 
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